Chapter 9 - ARIMA - Exercises

1. Figure 9.33 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

a. Explain the differences among these figures. Do they all indicate that the data are white noise?

The key difference between the three graphs is the blue lines, which represent \(\pm2 / \sqrt{T}\), where \(T\) is the length of the time series.

For each series, all of the points lie within these bounds, so they are all considered white noise.

b. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

The critical values are at different distances because they are inversely proportional to \(T\), the length of the time series.

The autocorrelations are different because, as \(T\) decreases and if the series is white noise, each lag will approach 0, the mean of the white noise.

2. A classic example of a non-stationary series is the daily closing IBM stock price series (data set fma::ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

fma::ibmclose %>% 
    as_tsibble() %>% 
    gg_tsdisplay(value, plot_type = 'partial')

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

When data have trends, the autocorrelations are positive and slowly decrease, which is what we see in the ACF. With the partial autocorrelation the linear dependence of previous lags removed, so we see that there is a strong positive relationship between the next day’s stock price and the previous day’s.

3. For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

a. `expsmooth::usnetelec`

uselec <-
  expsmooth::usnetelec %>%
  as_tsibble()

gg_tsdisplay(uselec, value, plot_type = 'partial')

There does not appear to be any seasonality with this series, so no Box-Cox transform is required. We perform a first difference:

uselec <-
  uselec %>% 
  mutate(value = difference(value))

gg_tsdisplay(uselec, value, plot_type = 'partial')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

The lags within the ACF and PACF plots are all less than the our white noise limits. We perform a unit root test, which is a hypothesis test. The null hypothesis is that the data is stationary.

uselec %>% 
  features(value, unitroot_kpss)

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.158         0.1

With a \(p\) value of 0.1 - and assuming a critical value of 0.05 - we can’t reject the null hypothesis, and we conclude the data is therefore stationary.

b. United states GDP from global_economy

us_economy <-
  global_economy %>%
  filter(Country == 'United States') 

autoplot(us_economy, GDP)

There is no variable seasonality, however the trend increases with time. A Box-Cox should help here:

us_economy %>% 
  features(GDP, guerrero)

## # A tibble: 1 x 2
##   Country       lambda_guerrero
##   <fct>                   <dbl>
## 1 United States           0.282

us_economy %>% 
  autoplot(box_cox(GDP, guerrero(GDP)))

us_economy %>% gg_tsdisplay(difference(box_cox(GDP, guerrero(GDP))), plot_type = 'partial')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

us_economy %>% gg_tsdisplay(box_cox(GDP, guerrero(GDP)) %>% difference() %>% difference(), plot_type = 'partial')

## Warning: Removed 2 row(s) containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_point).

The first difference does not create a statonary series, but the second difference appears to.

c. expsmooth::mcopper

mcopper <-
  expsmooth::mcopper %>% 
  as_tsibble()

gg_tsdisplay(mcopper, value, plot_type = 'partial')

This does not look seasonal, but does look cyclic with a trend upwards, hence non-stationary. We see the trend in the ACF decreasing exponentially. There is also some increased variation, so a Box-Cox transform may be beneficial:

bcg <- function(x) { box_cox(x, guerrero(x)) }

mcopper %>% autoplot(bcg(value))

mcopper %>% gg_tsdisplay(difference(bcg(value)), plot_type = 'partial')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

mcopper %>% gg_tsdisplay(difference(difference(bcg(value))), plot_type = 'partial')

## Warning: Removed 2 row(s) containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_point).

First and second difference look stationary, however there still appears to be some correlation (and negative correlation) in the first lags. We try a KPSS test:

mcopper %>% 
  features(difference(bcg(value)), unitroot_kpss)

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0464         0.1

At the first difference of the Box-Cox transform it appears we can’t reject the null hypothesis that the series is stationary.

d. expsmooth::enplanements

enp <-
  expsmooth::enplanements %>% as_tsibble()

autoplot(enp, value)

enp %>% gg_tsdisplay(bcg(value), plot_type = 'partial')

This series shows seasonality with increasing variance, and an upward trend.

First off would be a Box-Cox transformation and a first difference:

enp %>% gg_tsdisplay(difference(bcg(value)), plot_type = 'partial')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

The lags at six and 12 months are still visible. We’ll change the first difference to a 6 months:

enp %>%
  mutate(
    bc_value = bcg(value),
    bc_seasonal_diff = difference(bc_value, lag = 12),
    bc_second_diff = difference(bc_seasonal_diff)
  ) %>% 
  gg_tsdisplay(bc_second_diff, plot_type = 'partial')

## Warning: Removed 13 row(s) containing missing values (geom_path).

## Warning: Removed 13 rows containing missing values (geom_point).

e. expsmooth::visitors

visitors <- as_tsibble(expsmooth::visitors)

gg_tsdisplay(visitors, value, plot_type = 'partial')

Increasing variance suggests a Box-Cox, and the PACF suggests a difference of 1 and 12.

visitor_difference <-
  visitors %>% 
  mutate(
    value_bc = bcg(value),
    first_difference = difference(value_bc),
    second_difference = difference(first_difference, lag = 12)
  )

visitor_difference %>% gg_tsdisplay(second_difference)

## Warning: Removed 13 row(s) containing missing values (geom_path).

## Warning: Removed 13 rows containing missing values (geom_point).

## Warning: Removed 13 row(s) containing missing values (geom_path).

visitor_difference %>% features(second_difference, unitroot_kpss)

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0172         0.1

4. For the expsmooth::enplanements data, write down the differences you chose above using backshift operator notation.

We had a 12 month difference, followed by a one month difference. The 12 month difference is \(y_t - y_{t-12} = (1 - B^{12})y_t\), and the 1 month lag is \(y_t - y_{t-1} = (1 - B)y_t\). So we end up with \((1 - B^{12})(1 - B)y_t\).

5. For your retail data (from Exercise 6 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

set.seed(1232680)
myseries <- aus_retail %>%
 filter(`Series ID` == sample(aus_retail$`Series ID`,1))

autoplot(myseries, Turnover)

We’ve got increasing variability, so we’ll do a Box-Cox. The data is seasonal, so we’ll perform a 12 month difference, then a single lag difference.

myseries <-
  myseries %>% 
  mutate(Turnover_Diff = difference(difference(bcg(Turnover)), lag = 12))

gg_tsdisplay(myseries, Turnover_Diff, plot_type = 'partial')

## Warning: Removed 13 row(s) containing missing values (geom_path).

## Warning: Removed 13 rows containing missing values (geom_point).

6. Use R to simulate and plot some data from simple ARIMA models.

a. Use the following R code to generate data from an AR(1) model with $\phi_1=0.6$ and $\sigma^2=1$. The process starts with $y_1=0$

ar1 <- function(theta = 0.6) {
  y <- numeric(100)
  e <- rnorm(100)
  for(i in 2:100)
    y[i] <- theta*y[i-1] + e[i]
  
  a <- tsibble(idx = seq_len(100), y = y, index = idx)
  return(a)
}

b. Produce a time plot for the series. How does the plot change as you change $\phi_1$?

autoplot(ar1(0), y)

autoplot(ar1(), y)

autoplot(ar1(.90), y)

With a \(\phi = 0\), the process is purely random white noise. As \(\phi\) gets bigger, we see more influence of the previous values, and more wandering in a single direction for periods of time.

c. Write your own code to generate data from an MA(1) model with $\theta=0.6 and $\sigma^2=1$

ma1 <- function(theta, n = 100L) {
  y <- numeric(n)
  e <- rnorm(n)
  for(i in 2:n)
    y[i] <- theta * e[i-1] + e[i]
  
    tsibble(idx = seq_len(n), y = y, index = idx)
}

d. Produce a time plot for the series. How does the plot change as you change $\theta_1$?

ma1(0) %>% autoplot(y)

ma1(.5) %>% autoplot(y)

ma1(1) %>% autoplot(y)

ma1(10) %>% autoplot(y)

Changing the value of \(\theta\) doesn’t change the overall appearance of the plot, i.e. it still appears to be white noise. However it does increase the magnitude of the swings in each direction.

e. Generate data from an ARMA(1,1) model with $\phi _1=0.6$, $\theta_1=0.6$ and $\sigma^2=1$.

arma_1_1 <- function(phi, theta, n = 100L) {
  y <- numeric(n)
  e <- rnorm(n)
  
  for (i in 2:n) {
    y[i] <- (phi * y[i - 1]) + (theta * e[i - 1]) + e[i]
  }
  
  tsibble(idx = seq_len(n), y = y, index = idx)
}

arma <- arma_1_1(phi = .6, theta = .6)

f. Generate data from an AR(2) model with ϕ1=−0.8, ϕ2=0.3 and σ2=1. (Note that these parameters will give a non-stationary series.)

ar2 <- function(theta_1, theta_2, n = 100L) {
  y <- numeric(n)
  e <- rnorm(n)
  for(i in 3:n)
    y[i] <- theta_1*y[i-1] + theta_2*y[i-2] + e[i]
  
  tsibble(idx = seq_len(100), y = y, index = idx)
}

g. Graph the latter two series and compare them.

arma_1_1(.6, .6) %>% autoplot(y)

ar2(-.8, .3) %>% autoplot(y)

The AR2 model has an oscillation in it, so it is not stationary.

7. Consider fpp2::wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.

wmurders <-
  fpp2::wmurders %>% 
  as_tsibble()

a. By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data.

wmurders %>% autoplot(value)

The first thing to notice is that we don’t have a stationary time series. There is no seasonality, but it does wander up and down for significant periods. We perform difference, and unit test to determine if the series is stationary.

wmurders <-
  wmurders %>% 
  mutate(first_diff = difference(value))

gg_tsdisplay(wmurders, first_diff, plot_type = 'partial')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

wmurders %>% features(first_diff, unitroot_kpss)

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.470      0.0485

The series looks reasonably stationary, and our prominent lags are only just outside the critical values. We’re just under the level to reject the null hypothesis, but we’ll stick with \(d = 1\) in the ARIMA model.

Looking at the lags, lag 2 has significance in both the ACF and PACF plot. We’ll start with an \(ARIMA(2,1,0)\)

a. Should you include a constant in the model? Explain    .

From Rob’s Blog:

the inclusion of a constant in a non-stationary ARIMA model is equivalent to inducing a polynomial trend of order \(d in\) the forecast function. (If the constant is omitted, the forecast function includes a polynomial trend of order \(d−1\).) When \(d=0\), we have the special case that \(\mu\) is the mean of \(y_t\).

Looking at the first difference, I beleive there is a linear trend downwards. Adding in a constant intrduces a trend of order 1 (linear), thus addition of a constant is warranted.

Correct answer: A constant would imply a drift in the original data which does not look correct, so we omit the constant.

a. Write this model in terms of the backshift operator.

$$ (1 - B)y_t = c + _1 B(1 - B)y_t + _2 B^2(1 - B)y_t + _t \

(1 = B)y_t = c + _1 (B - B^2)y_t + _2 (B^2 - B^3)y_t + _t $$

a. Fit the model using R and examine the residuals. Is the model satisfactory?

wm_arima <-
  wmurders %>% 
  model(ARIMA(value ~ pdq(2,1,0)))

report(wm_arima)

## Series: value 
## Model: ARIMA(2,1,0) 
## 
## Coefficients:
##           ar1     ar2
##       -0.0572  0.2967
## s.e.   0.1277  0.1275
## 
## sigma^2 estimated as 0.04265:  log likelihood=9.48
## AIC=-12.96   AICc=-12.48   BIC=-6.99

gg_tsresiduals(wm_arima)

- The lag plot shows all lags within the critical value. - No patterns in the residuals, some increasing variance. - Residuals are appoximately normally distributed.

a. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.

First we programmatically forecast:

wm_arima %>% 
  forecast(h = 3) %>% 
  pull(.mean)

## [1] 2.553355 2.533802 2.524231

We have an ARIMA(2,1,0), so our our forecasting equation will be:

\[ y^\prime_t = \phi_1 By^\prime_t + \phi_2 B^2 y^\prime_t \] As we have differenced once, we have \(y^\prime_t = (1 - B)y_t\). We substitute this in:

$$ (1 - B)y_t = _1 B(1-B)y_t + _2 B^2(1 - B)y_t + _t\

(1 - B)y_t = _1 (B - B^2)y_t + _2 (B^2 - B^3)y_t + _t \

y_t - y_{t-1} = 1(y{t-1} - y_{t-2}) + 2(y{t-2} - y_{t-3}) + _t \

y_t = y_{t-1} + 1(y{t-1} - y_{t-2}) + 2(y{t-2} - y_{t-3}) + _t $$

phi <- wm_arima %>% tidy() %>% pull(estimate)
# Add one so we can do 'n - 1' to be the same as our 't - 1' in the equation.
n <- nrow(wmurders) + 1
y <- wmurders[['value']]

for (x in 1:3) {
  n <- length(y) + 1
  y[n] <- y[n - 1] + (phi[1] * (y[n-1] - y[n-2])) + (phi[2] * (y[n-2] - y[n-3]))
  print(y[n])
}

## [1] 2.553355
## [1] 2.533802
## [1] 2.524231

We can see that this matches the output from our forecast() function.

a. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.

wm_arima %>% 
  forecast(h = 3) %>% 
  autoplot(slice_tail(wmurders, n = 20))

a. Does ARIMA() give the same model you have chosen? If not, which model do you think is better?

wmurders %>% 
  model(ARIMA(value)) %>% 
  report()

## Series: value 
## Model: ARIMA(1,2,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.2434  -0.8261
## s.e.   0.1553   0.1143
## 
## sigma^2 estimated as 0.04632:  log likelihood=6.44
## AIC=-6.88   AICc=-6.39   BIC=-0.97

The ARIMA() function gives us a \(ARIMA(1,2,1)\) model. This has a much lower (and therefore much better AIC that the \(ARIMA(2,1,0)\) that I chose.

8. Consider fpp2::austa, the total international visitors to Australia (in millions) for the period 1980-2015.

a. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.

austa <-
  fpp2::austa %>% as_tsibble()

autoplot(austa, value)

austa_arima <-
  austa %>% 
  model(ARIMA(value)) 

austa_arima %>% report()

## Series: value 
## Model: ARIMA(0,1,1) w/ drift 
## 
## Coefficients:
##          ma1  constant
##       0.3006    0.1735
## s.e.  0.1647    0.0390
## 
## sigma^2 estimated as 0.03376:  log likelihood=10.62
## AIC=-15.24   AICc=-14.46   BIC=-10.57

The ARIMA() function has chosen an \(ARIMA(0,1,1)\)

Looking at the residuals:

austa_arima %>% gg_tsresiduals()

The residuals look like white noise, however they don’t look as normal as I’d like.

Forcasting:

austa_arima %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

a. Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part a. Remove the MA term and plot again.

# ARIMA(0,1,1) with no drift
austa %>% 
  model(ARIMA(value ~ 0 + pdq(0,1,1))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

# ARIMA(0,1,0) with no drift
austa %>% 
  model(ARIMA(value ~ 0 + pdq(0,1,0))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

Removing the dift in the differenced series removes the trend in the original, which has a very large effect on the forecasts.

a. Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.

# ARIMA(2,1,3) with drift
austa %>% 
  model(ARIMA(value ~ 1 + pdq(0,1,0))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

# Constant removed
austa %>% 
  model(ARIMA(value ~ 0 + pdq(0,1,0))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

Again, without the drift, the trend is not seen in the forecasts.

a. Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.

# ARIMA(0,0,1) with drift
austa %>% 
  model(ARIMA(value ~ 1 + pdq(0,0,1))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

# Constant removed
austa %>% 
  model(ARIMA(value ~ 1 + pdq(0,0,0))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

A stationary model with a constant has long term forecasts equal to the mean.

a. Plot forecasts from an ARIMA(0,2,1) model with no constant.

# Constant removed
austa %>% 
  model(ARIMA(value ~ 0 + pdq(0,2,1))) %>% 
  forecast(h = 10) %>% 
  autoplot(austa)

This model visually appears to do well and captures the trend, even though there isn’t a constant. From the notes we recall that when \(c = 0\) and \(d = 2\), the long term forecasts will follow a straight line.

9. For the United States GDP series (from global_economy):

   1. if necessary, find a suitable Box-Cox transformation for the data;

   2. fit a suitable ARIMA model to the transformed data using ARIMA();

   3. try some other plausible models by experimenting with the orders chosen;

   4. choose what you think is the best model and check the residual diagnostics;

   5. produce forecasts of your fitted model. Do the forecasts look reasonable?

   6. compare the results with what you would obtain using ETS() (with no transformation).

10. Consider fpp2::austourists, the quarterly number of international tourists to Australia for the period 1999–2010.

    1. Describe the time plot.
    2. What can you learn from the ACF graph?
    3. What can you learn from the PACF graph?
    4. Produce plots of the seasonally differenced data (1−B4)Yt
    5. What model do these graphs suggest?
    6. Does ARIMA() give the same model that you chose? If not, which model do you think is better?
    7. Write the model in terms of the backshift operator, then without using the backshift operator.

11. Consider fpp2::usmelec, the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 – June 2013). In general there are two peaks per year: in mid-summer and mid-winter.

    1. Examine the 12-month moving average of this series to see what kind of trend is involved.
    2. Do the data need transforming? If so, find a suitable transformation.
    3. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
    4. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
    5. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
    6. Forecast the next 15 years of electricity generation by the U.S. electric industry. Get the latest figures from the EIA to check the accuracy of your forecasts.
    7. Eventually, the prediction intervals are so wide that the forecasts are not particularly useful. How many years of forecasts do you think are sufficiently accurate to be usable?

12. For the expsmooth::mcopper data:

    1. if necessary, find a suitable Box-Cox transformation for the data;
    2. fit a suitable ARIMA model to the transformed data using ARIMA();
    3. try some other plausible models by experimenting with the orders chosen;
    4. choose what you think is the best model and check the residual diagnostics;
    5. produce forecasts of your fitted model. Do the forecasts look reasonable?
    6. compare the results with what you would obtain using ETS() (with no transformation).

13. Choose one of the following seasonal time series: the Australian production of electricity, cement, or gas (from aus_production).

    1. Do the data need transforming? If so, find a suitable transformation.
    2. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
    3. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
    4. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
    5. Forecast the next 24 months of data using your preferred model.
    6. Compare the forecasts obtained using ETS().

14. For the same time series you used in the previous exercise, try using a non-seasonal model applied to the seasonally adjusted data obtained from STL. Compare the forecasts with those obtained in the previous exercise. Which do you think is the best approach?

15. For the Australian tourism data (from tourism):

    1. Fit a suitable ARIMA model for all data.
    2. Produce forecasts of your fitted models.
    3. Check the forecasts for the “Snowy Mountains” and “Melbourne” regions. Do they look reasonable?

16. For your retail time series (Exercise 5 above):

    1. develop an appropriate seasonal ARIMA model;
    2. compare the forecasts with those you obtained in earlier chapters;
    3. Obtain up-to-date retail data from the ABS website (Cat 8501.0, Table 11), and compare your forecasts with the actual numbers. How good were the forecasts from the various models?

17. Consider fma::sheep, the sheep population of England and Wales from 1867–1939.

    1. Produce a time plot of the time series.

    2. Assume you decide to fit the following model: yt=yt−1+ϕ1(yt−1−yt−2)+ϕ2(yt−2−yt−3)+ϕ3(yt−3−yt−4)+εt, where εt is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?

    3. By examining the ACF and PACF of the differenced data, explain why this model is appropriate.

    4. The last five values of the series are given below:

    Year 1935 1936 1937 1938 1939 Millions of sheep 1648 1665 1627 1791 1797

    The e stimated parameters are ϕ1=0.42 ϕ2=−0.20, and ϕ3=−0.30

    Without using the forecast function, calculate forecasts for the next three years (1940–1942).

    5. Now fit the model in R and obtain the forecasts using forecast. How are they different from yours? Why?

18. The annual bituminous coal production in the United States from 1920 to 1968 is in data set fma::bicoal.

    1. Produce a time plot of the data.

    2. You decide to fit the following model to the series: yt=c+ϕ1yt−1+ϕ2yt−2+ϕ3yt−3+ϕ4yt−4+εt where yt is the coal production in year t and εt is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?

    3. Explain why this model was chosen using the ACF and PACF.

    4. The last five values of the series are given below.

    Year 1964 1965 1966 1967 1968 Millions of tons 467 512 534 552 545

    The estimated parameters are c=162.00, ϕ1=0.83, ϕ2=−0.34, ϕ3=0.55, and ϕ4=−0.38. Without using the forecast function, calculate forecasts for the next three years (1969–1971).

5. Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?

19. Before doing this exercise, you will need to install the Quandl package in R.

    1. Select a time series from Quandl. Then copy its short URL and import the data using y <- as_tsibble(Quandl("?????", api_key="?????"), index = Date)

    2. Plot graphs of the data, and try to identify an appropriate ARIMA model.

    3. Do residual diagnostic checking of your ARIMA model. Are the residuals white noise?

    4. Use your chosen ARIMA model to forecast the next four years.

    5. Now try to identify an appropriate ETS model.

    6. Do residual diagnostic checking of your ETS model. Are the residuals white noise?

    7. Use your chosen ETS model to forecast the next four years.

    8. Which of the two models do you prefer?